Stability of negative ionization fronts: regularization by electric screening? 
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We recently have proposed that a reduced interfacial model for streamer propagation is able 
to explain spontaneous branching. Such models require regularization. In the present paper we 
investigate how transversal Fourier modes of a planar ionization front are regularized by the electric 
screening length. For a fixed value of the electric field ahead of the front we calculate the dispersion 
relation numerically. These results guide the derivation of analytical asymptotes for arbitrary fields: 
for small wave- vector k, the growth rate s(k) grows linearly with k, for large k, it saturates at some 
positive plateau value. We give a physical interpretation of these results. 



I. INTRODUCTION 

Streamers generically appear in electric breakdown 
when a sufficiently high voltage is suddenly applied to 
a medium with low or vanishing conductivity. They con- 
sist of extending fingers of ionized matter and are ubiq- 
uitous in nature and technology. Frequently they are 
observed to branch 0, 0] • There is a traditional qualita- 
tive concept for streamer branching based on rare photo- 
ionization events 0, El IS IE 0- However, our recent work 
(E IE E] nas shown that even the simplest, fully de- 
terministic streamer model without photo-ionization can 
exhibit branching. In particular, we have proposed |8| 
that a streamer approaching the Lozansky-Firsov limit 
of ideal conductivity |Tlj can branch spontaneously due 
to a Laplacian interfacial instability 12]. This mecha- 
nism is quite different from the one proposed previously. 
It requires less microscopic physical interaction mecha- 
nisms, but is based on a dynamically evolving internal 
interfacial structure of the propagating streamer head. 
Analytical branching predictions from the simplest type 
of interfacial approximation can be found in |l(|. 

However, the simple interfacial model investigated in 
|lfj| requires regularization to prevent the formation of 
cusps. The nature of this regularization has to be de- 
rived from the underlying gas discharge physics; it re- 
cently has been subject of debate [l3l,ll4j . We argue that 
one regularization mechanism is generically inherent in 
any discharge model, namely the thickness of the electric 
screening layer. This is the subject of the present pa- 
per: we study how the electric screening layer present in 
the partial differential equations of the electric discharge 
influences the stability of an ionization front, correc ting 
the simple interfacial model proposed in [E Ej Ei Il5j 
and solved in El- To be precise, we derive the disper- 
sion relation for transversal Fourier-modes of a planar 
ionization front. We treat a negative front in a model as 

in H HE EE [13], but with vanishm S electron dif- 
fusion and under the assumption that the state ahead 
of the ionization front is completely non-ionized. We 
have sho wn p reviously that the analysis of the full model 
[EIEEjEjEjE] is mathematically nonstandard and 
challenging due to the "pulled" nature [l8L E] of the 



front. Pulling is a mode of front propagation where the 
spatially half- infinite leading edge of a front dominates its 
behavior. However, for vanishing electron diffusion and 
propagation into a non-ionized state, the leading edge of 
the ionization front is completely eliminated and replaced 
by a discontinuous jump of the electron density to some 
finite value. This corresponds to the fact that neglecting 
electron diffusion changes the equation of electron mo- 
tion from parabolic to hyperbolic type. Putting D e = 
in the present paper and considering propagation into a 
non-ionized state, we get rid of leading edge and pulling, 
but in turn we have to analyze discontinuous fronts. 

Here we anticipate the result of the paper: if the field 
far ahead of a planar negative ionization front is Eoo, 
then a transversal Fourier perturbation with wave vector 
k grows with rate 
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where a(E) is the effective impact ionization coefficient 
within a local field E; a sets the size of the inverse electric 
screening length. The behavior for large A; is a correction 
to the interfacial model treated in [l(| ; in that model we 
would have s(k) = \Eoo\k for all k. The asymptotes Q 
have been quoted already in [E Il5| . however, without 
derivation. Their derivation based on numerical results 
and asymptotic analysis together with a discussion of the 
underlying physical mechanisms are the content of the 
present paper. 

In detail, the paper is organized as follows: in Sec. II 
we summarize the minimal streamer model in the limit of 
vanishing diffusion and recall multiplicity, selection and 
analytical form of uniformly translating planar front so- 
lutions; we then derive the asymptotic behavior at the 
position of the shock and far behind the shock, and we 
discuss two degeneracies of the problem. In Sec. Ill we 
set up the framework of the linear perturbation analysis 
for transversal Fourier modes, first the equation of mo- 
tion and then the boundary conditions and the solution 
strategy. In Sect. IV we present numerical results for the 
dispersion relation for field E^ = — 1, and we derive the 
asymptotes Q analytically for arbitrary E^ . The small 
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k limit is related to one of the degeneracies of the unper- 
turbed problem, for the large k limit we also present a 
physical interpretation. Sect. V contains conclusions and 
outlook. 



II. MINIMAL STREAMER MODEL AND 
PLANAR FRONT SOLUTIONS 

A. The minimal model 

We investigate the minimal streamer model, i.e., a 
"fluid approximation" with local field-dependent impact 
ionization reaction in a non-attaching gas like argon or 
nitrogen 0, H [H [H GJ, El HI For physical param- 
eters and dimensional analysis, we refer to our previous 
discussions in [H, 0, 0, 0] . When electron diffusion is 
neglected (D e = 0) , the dimensionless model has the form 



d t a - V • (a E) = a /(E) , 
d tP = a /(E) , 
V-E = p-o , E 



-v</>, 



(2) 
(3) 
(4) 



where a is the electron and p the ion density and E the 
electric field. Here the electron current is assumed to 
be crE, and the ion current is neglected. Electron-ion 
pairs are assumed to be generated with rate er/(E) = 
a|E| a(|E|) where tr|E| is the absolute value of electron 
current and a(E) the effective impact ionization cross 
section within a field E. Hence /(E) is 



/(E) = |E| a(|E| 



(5) 



For numerical calculations, we use the Townsend approx- 
imation 



a(|E|) 



(6) 



For analytical calculations, an arbitrary function a(E) 
can be chosen where we only assume that 



/(E) = /(|E|) and a(0) = 



(7) 



= /'(0). For 
does not decrease 



The last identity entails that /(0) 
certain results we also need that a(|E 
when |E| increases, hence that a' > 0. 

Note that the electrons are the only mobile species 
and the source of additional ionization, while ion density 
p and electric potential <f> or field E follow the dynamics 
of the electron density er, and couple back onto it. 



B. Uniformly translating ionization fronts: 
analytical solutions and multiplicity 

We now recall essential properties of uniformly trans- 
lating planar front solutions of Eqs. @~© and {7|). First 
of all, a constant mode of propagation requires a planar 
density distribution that we assume to vary only in the z 



direction: (er, p) = (a(z,t), p(z,t)); the particle densities 
for large positive z are assumed to vanish. The field far 
ahead of the front in the non-ionized region at z — » oo 
has to be constant in time and as a consequence of Q 
also constant in space: 



E 







+oo 

— oo 



(8) 



where z is the unit vector in z direction. For the bound- 
ary condition at z — > — oo we assumed that the ionized 
region behind the front extends to — oo. This implies, 
that a fixed amount of charge J(p — cr)dz = Eqo is trav- 
eling within the front according to {3J and ©, and no 
currents flow far behind the front in the ionized and elec- 
trically screened region. 

For the further analysis, a coordinate system (x, y, £ = 
z — vt) moving with velocity v in the z direction is used. 
Then the equations |[5J Q read 



<9t<7 — vd^a 
d t p - vd^p 



(p - a) a + (Va) • (V0) - a/(|V0|) = , 
<t/(|W|)=0, 



p-a + V 2 (f> = , 



(9) 



where we expressed all quantities by electron density a, 
ion density p and electric potential 4>- 

A front propagating uniformly with velocity v is a so- 
lution of JSJl, @ where er, p and </> depend of £ only. With 

Vc6 = d^cf) z — —E z, such a front solves 

{v + E)d 6 a + (p - a) a + af(\E\) = 0, (10) 
vdtp + <rf(\E\) = , (11) 
p-o-d^E = . (12) 

For use in the later sections, we now briefly recall the 
analytical solutions [|2j of these equations. Subtract pif) 
from UlUfl . use l|12l) to get a complete differential, inte- 
grate and use l|12|) again to get —vd^E + aE = const.. 
The integration constant is fixed by the condition E — > 
at £ — > — oo from Eq. JSJl, and we find 



idt-E + aE = 



(13) 



The front equations then reduce to two ordinary differ- 
ential equations for er and E 

dt[(v + E)a] = -af(E) , f(E) = \E\a(E), 
vd s ln\E\ = er, (14) 

that can be solved analytically: 

°\ E \ = ^T-5P[E\, (15) 
P\E] 



v + E 
' B ~'/(*) 



dx 



6-6 



E\ x J\E 

E{ ^ 2) v + x dx 



a(x)dx, (16) 
(17) 
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This gives us a and p as functions of E, and the space 
dependence E — E(£) implicitly as £ = £(E) in the last 
equation. It follows immediately from Ijl7|l that E(£) is a 
monotonic function, and hence that the space charge q — 
p — a = d^E has the same sign everywhere. According 
to Ijl6|l . p(£) is a monotonic function, too. 

Up to now, the front velocity v as a function of the 
asymptotic field Eoo is not yet fixed. Indeed for any non- 
vanishing far field E^, there is a continuous family of 
uniformly translating front solutions parametrized by v 
E2- HII > since the front propagates into an unstable state 
jl8j. In particular, for E^ > there is a dynamically sta- 
ble solution for any velocity v > 0, and for E^ < 0, there 
is a dynamically stable solution for any v > \Eoo\- These 
bounds on v can be derived directly from Eqs. (|15|l l|17|) 
with boundary condition (jSJ) and the condition that the 
densities a and p are non-negative for all £. 

This continuous family of solutions parametrized by 
v is associated with an exponentially decaying electron 
density profile in the leading edge 0, 0] : an electron 
profile that asymptotically for large £ decays like <r(£) oc 
e~ A £ with A > 0, will propagate with velocity 



-E n 



f(Eoo) 

A 



in a field E^ < 0. 



(18) 



It will "pull" an ionization front along with the same 
speed. (For E^ > 0, the same equation applies for all 
A > /(-Boo) /Eoo, hence for v > 0). 



C. Dynamical selection of the shock front solution 
and its particular properties 

In practice, not all these uniformly propagating solu- 
tions are observed as asymptotic solutions of the full dy- 
namical problem but only a specific one that is 
called the selected front. For a negative ionization front, 
it propagates with the velocity [12j 



v= \E B 



for Eoa < 



(19) 



The selection takes place through the initial conditions 
[l8[ : If the electron density strictly vanishes beyond a 
certain point £0 at time t = 



a = = p for £ > £ at t = 0, 



(20) 



then this stays true for all later times t in the comov- 
ing frame £. Only initial conditions that decay expo- 
nentially like e -A ^ for £ — > 00, approach a solution with 
the larger velocity (|18|) . Such an exponential decay is a 
very specific initial condition; furthermore, such a lead- 
ing edge will generically be cut off for very small densities 
by the physical break-down of the continuum approxima- 
tion. Therefore the physically relevant solution is the one 
with velocity i(T§|) and absent leading edge as in (|2*U|l . 
The complete absence of the leading edge (A = 00) is 
generic for the hyperbolic equation J2J), i.e., for vanish- 
ing electron diffusion. We will restrict the analysis of 



fronts and their linear perturbations to propagation into 
a completely non-ionized state (|2U[) in the remainder of 
the paper. 

In contrast to all other uniformly translating fronts 
with v > —Eqo, the selected front with v = —Eoq exhibits 
a discontinuity of the electron density at some point £ 
which corresponds to v + E(£) — > 0. We choose the co- 
ordinates such that the discontinuity is located at £ = 0. 
The situation is shown in Fig. ^ for a uniformly translat- 
ing front with velocity v = 1 within a far field E^ = — 1 . 




FIG. 1: Electron density a (solid line in first plot), ion density 
p (dotted line in first plot) and electric field E (second plot) 
for a negative ionization shock front moving with v = \Eoo\ 
in the comoving frame £ = z — vt. The far field is = — 1. 



A discontinuity of a means that d^a is singular at this 
position. On the other hand, the expression a(p — a + 
f(E)) in Eq. l|lUfl is finite or vanishing, therefore the 
product (v + E)d^a in Eq. (|l(Jfl may not diverge either. 
Hence (v + E) has to vanish at the position of the discon- 
tinuity, and therefore E — E^ = —v at the position of 
the front. Furthermore, since (v + E) — > for £ f while 
d^a is bounded for £ < — as we will derive explicitly 
below in Eq. I|28() — we have 



lim 



e(o a c( r = o. 



(21) 



The fact that c(£) in Fig. increases monotonically 
up to the position of the shock, is generic and can be 
seen as follows: according to (I1U|I . and since (v + E) > 
and a > 0, the sign of d^a is identical to the sign of 
a — p — f(E). With the help of the exact solutions l|15|) 
and l|lbfl . with the definition of f(E) in JSJ and with 
identifying v = [E^L we find 



a-p-f(E) = \E\ 



\E\ 



a(x) — a(E) 
-\E\ 



dx > 0. 



(22) 



So er(£) increases monotonically for growing £ up to £ = 
as long as a{E) increases monotonically with E. This is 
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the case for Townsend form © or more generally for any 
a(E) that is monotonically increasing with E. 



D. Asymptotics near the shock front 

We now derive explicit expressions for er(£) etc. near 
the discontinuity. On approaching the position of the 
ionization shock front from below f f 0, the quantity 



e = v + E= \E a 



\E\ 



(23) 



is a small parameter. The ion density (|16[1 at this point 
can be expanded as 



p[E] = p[v - e] = a(v)e - a'(v) - + 0(e 3 ). (24) 

As the electron density is related to the ion density 
through cr[E] — p[E] v/e according to 115|) . it is 



a[E] = va(v) - va'(v)~ + 0(e 2 ). 

Eq. JTJJ evaluated for E(^ 2 = 0) = Eoo < reads 
v — x dx f e y dy 



(25) 



E(®\ PN x J p[v-y] v-y 



, (26) 



where in the last expression, the parameter e (I23|l is in- 
troduced. Insertion of l|24(l now yields an explicit relation 
between £ and E: 



■t = 



+ 0(e 2 ) 



va(v) 

or t = -va(v)£ + 0(£ 2 ). 



(27) 



Insertion of this approximately linear relation between 
e and £ into (|24|l and (|25|l together with the notation 
f(v) = va(v) results in 



= /(«) + 



^^C + 0(a)(28) 



P(0 = ( - f(v) a(v) £ + 0(e)) , (29) 

-25(0 - v + 6(-0(f(v)£ + O(e)), (30) 

where we used v — |i?oo| and the step function 0(x), 
defined as 0(x) = 1 or for x > or x < 0, respectively. 

E. Asymptotics far behind the shock front 

Far behind the front in the ionized region £ — > — oo, the 
fields approach lim^_ i _ 00 (a, p, E) = (er~ , p~ , E~) with 



= p~ = \ a(x) dx , E~ = . (31) 
Jo 



Expanding about this point as cr(£) = cr + cri(£) etc., 
we derive in linear approximation 



(32) 





with A given by 



x= a - = 



, . dx 
a(x) — . 



(33) 



Two eigenvalues of the matrix in (|32ll vanish. The third 
eigenvalue of the matrix is the positive parameter A, it 
produces the eigendirection 




for £ — » — oo, 



(34) 



that describes the asymptotic solution deep in the ionized 
region. The free parameter A > accounts for transla- 
tion invariance. 



F. Two degeneracies of the shock front 

We have fixed the initial condition (12011 and hence we 
have selected the front speed v = —E^. Therefore the 
degeneracy of solutions related to the profile of the lead- 
ing edge is removed. Still there are two degeneracies re- 
maining in the problem. The first one is the well known 
mode of infinitesimal translation that corresponds to the 
arbitrary position of the front. The second one is spe- 
cific for the present problem and will play a central role 
in the derivation of the analytical asymptote for small k 
in Sect. IV. It is the mode of infinitesimal change of far 
field E^. It corresponds to the arbitrariness of the field 
Eqo in the non-ionized region with a = = p ahead of 
the front and to the corresponding arbitrariness of the 
asymptotic ionization level a = a~ = p behind the front 
where the field vanishes. To set the stage for the later 
analysis, the necessary properties of the modes are given. 

An infinitesimal translation of the front in space gen- 
erates the linear mode (cr u pt, E t ) = (d^a, d^p, d^E) 



{v + E)d s a t = [2a- p- fja t -ap t + [af' -d e ajE t 
vd^pt = —fat+af'Et 

d € E t = pt-at (35) 

with the definition /' = d x f(\x\), so that f(E + E t ) — 
f - f'E t + ... for 15 < 0. With the notation ip t = —E t , 
the equations can be written in matrix form as 




No(0 




(36) 



5 



/ 2a - f 



N (0 



E 



v + E 



V 



v 
1 



d i a-af\ 
v + E 

-of 







(37) 



Note that the matrix Note) reduces to the matrix in Eq. 
P2| for f -» -oo, since (2a-f-p)/(v+E) -> o-/« = A 
etc. The limiting value of the vector (a t , p t ,tp t ) for £ — > 
is according to ffi fy -jlfl ty 



(38) 



The second mode is generated by an infinitesimal 
change of the far field Eoo and consecutively by an in- 
finitesimal change of the velocity v. The discontinuity is 
taken at the position £ — 0. In linear order, this variation 
creates a mode 




(T£ f = hm etc. 

c^o e 

that solves the inhomogeneous equation 



(39) 



No(£) 





dp /(v + E) 
d^p/v 


' (40) 

The inhomogeneity vanishes at £ — > — oo. Hence 
like the front solution itself and like the infinitesimal 
translation mode, also this mode has the eigendirec- 
tion (Sa~ , 5a~ , 0) + A (A, 0, 1) e A? + . . . asymptotically 
for £ — > — oo. The value of 5cr _ is given by 6a~ = 
da~/d\Eoo\ = a(E QO ) according to (J3TJ. For £ f 0, the 
limiting values of the fields are 



(41) 



which is the derivative of i|28|) - (|3U[) with respect to v 
l^oolate^O. 



III. SET-UP OF LINEAR STABILITY 
ANALYSIS 

We now can proceed to study the stability of a planar 
ionization shock front. The front propagates into the z 
direction. The perturbations have an arbitrary depen- 
dence on the transversal coordinates x and y. Within 
linear perturbation theory, they can be decomposed into 
Fourier modes. Therefore we need the growth rate s(k) of 
an arbitrary transversal Fourier mode to predict the evo- 
lution of an arbitrary perturbation. Because of isotropy 
within the transversal (x, y)-plane, we can restrict the 



analysis to Fourier modes in the x direction, so we study 
linear perturbations oc e st +' lkx . (The notation anticipates 
the exponential temporal growth of such modes.) 

In general, there can be a degeneracy of the dispersion 
relation s(k) for various profiles of the leading edge just 
as it is found also for the uniformly translating solutions 
in Section II. B. The constraint of a non-ionized initial 
condition (|20|) again will remove this degeneracy and fix 
s(k). In the present section, we will derive the equations 
and the boundary conditions for the Fourier modes. In 
Sect. IV, we will solve them numerically and derive the 
analytical asymptotes QJ. 



A. Equation of motion 

The linear perturbation theory could be set up within 
the coordinate system (x, £ = z — vt) that moves with the 
unperturbed constant velocity v = |i?oo|- This would, of 
course, lead to a set of equations that are linear in the 
perturbation. 

However, when the perturbation of a planar front 
grows, the position of the actual discontinuity of the elec- 
tron density will deviate from the position of the discon- 
tinuity of the unperturbed front. Within the coordinate 
system (x, £), this would lead to finite deviations within 
infinitesimal spatial intervals instead of infinitesimal de- 
viations within finite intervals. This conceptual difficulty 
can be avoided by formulating the perturbation theory 
within the coordinate system of the position of the per- 
turbed shock front (x, C) with 



C = £ - A k , £ = z-vt , A fe = 5 e 



ikx-\-st 



(42) 



where z is the rest frame, £ is the frame moving with 
the planar front and £ = marks the line of electron 
discontinuity of the actual front. Therefore we write the 
perturbation as 



a(x,C,t) = <7 (0 + ffi(0 A fc (x,t) 
p(x,C,t) = po(C) + Pi(0 A fc (x,t), 

0(a;,c,t) = MQ + MQ A k(x,t), 



(43) 



where ctq, po and </>o are the electron density, ion density 
and electric potential of the planar ionization shock front 
from the previous section. But these planar solutions 
here are shifted to the position of the perturbed front £. 
Therefore they do not move with their proper velocity 
v = —dt£, but with a slightly different velocity — dt( = 
v — sAfc. The price to pay is that the equations of the 
perturbation analysis become inhomogeneous, actually in 
a similar way as in |18| . The gain is that the derivation 
of the boundary conditions at the shock front becomes 
more comprehensible, and that later in Section V.B the 
identification of the analytical solution for small k with 
the mode (<te, pe^e) from the previous section becomes 
quite obvious. 
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Substitution of the expressions into © gives to 
leading order in the small parameter S 

(v + E ) <%<ti = (s + 2(7o — Po — f) &i 

-cr Pi + (d^aa - (J /') d ( (f>i - sd ( a , 
v d^pi = -f at + s pi - <7 f d C (f>i - sd ( p , 
(d 2 - k 2 ) 4> x = ax - p x + k 2 E . (44) 



Furthermore 



Here / = f(E ), f = d\ E \f(\E 



Eo 



and E = -c\0 o (O 

is the electric field of the uniformly translating front. As 
explained above, these equations are not completely lin- 
ear in (i7i, px,4>x), but contain the inhomogeneities sd^ao, 
sd^po and k 2 Eo. 

To elucidate the structure of Eq. I|44|l . we drop all in- 
dices and introduce the matrix notation 




M 



s,k 




sd c a/(v + E) 
sd^p/v 
-Ek 2 




/ s + 2a — f — p —a d^a — af 



M fl , fc (C) 



(45) 



V 



Here we introduced the auxiliary field 

ipx = d C (j)x, 



+ E 


v + E 


v + E 






s 


-of 





V 


V 


V 


1 


-1 





k 2 








1 


J 



(46) 



(47) 



that corresponds to the perturbation E\ of the electric 
field, but with reversed sign. 



B. Boundary conditions at the discontinuity 

Having obtained the perturbation equations, we are 
now in the position to derive the boundary conditions. 
First we consider the boundary conditions at £ = where 
we make explicit use of the initial condition JSUJ. The 
boundary conditions arise from the boundedness of the 
electron density to the left of the shock front at £ \ 0, and 
from the continuity of all other fields across the position 
C = of the shock front. 

As discussed in Section II. D, for the uniformly propa- 
gating shock front, the quantity (v + E) 8q<j vanishes as 
C T 0, since (v + E) vanishes and d^cr is bounded. Since 
this should hold both for the full solution as well as for 
the unperturbed solution, it also holds for the perturba- 
tion 



lim 

CTO 



v + E(C) d c ax = 0. 



(48) 



E(0 



drax =0 for ( > 0. 



(49) 



This identity is trivial for £ > 0, but nontrivial for £ = 0. 
When the explicit expressions (|28[1 - H30[) are inserted into 
(PHjl . we find 

(v + E)d c ax = {s + f{v))<Jx-f(v)px-f{v)f'{v)^x 
+(i>x-s)d ( a + O(0- (50) 

First of all, d^a is singular at £ = 0, since d^cr oc 
d(6(-() = -5(0- Therefore gHJ requires that the coef- 
ficient of d^a must vanish 



(51) 



which gives the first boundary condition. Second, apply- 
ing now Ij48(l yields the second boundary condition 

(s + f(v))ax(0) - f(v) Pl (0) - /(«)/'(t#i(°) = 0- ( 52 ) 

Due to the discontinuity, actually two boundary condi- 
tions JHU and (JSJJ result from (gEl and 

In a second step the continuity of the other fields across 
C = is evaluated. The continuity of p we get from Ijllf) 
and the fact, that a and |E| are bounded for all £. It 
immediately yields the third boundary condition 



Pi(0) = 0, 



(53) 



just like for the unperturbed equation. Finally, for the 
boundary conditions on field and potential, it is helpful 
that there is an exact solution for the non-ionized region 
at £ > for a boundary with the harmonic form II42II . 
Since ahead of the front there are no particles a = = p, 
there are also no space charges, and for the potential, 
one has to solve V 2 = with the limit E = — V</> — ► 
Eqo £ = — v £ as C —* oo. The general solution for £ > 
is 

4> = vZ + Sce~ ki e ikx+st 

= v( + 5{v + c e- kC ) e lkx+st + 0(5 2 ) (54) 

with the yet undetermined integration constant c. Here 
we chose the gauge 4>o(£, = 0) = for the unperturbed 
electric potential. 

Now <j) always is continuous, and E = — V</> is con- 
tinuous, because the charge density \p — cr| < oo in lf4"| 
everywhere. The continuity of <p at £ = implies 



</>i (0) = v + c, 



(55) 



the continuity of d x <p yields the same condition, and the 
continuity of d^4> implies 



V>i(0) = -ck. 



(56) 



The five boundary conditions {51 )1 -1(53 )1 and l)55 )l -l(5o )l 

determine the value of the integration constant 



s 
k 



(57) 
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in l|54|) and the values of the four fields at £ = 



//»«/(«)/(« + /(«)) 
Pi CTO 

V (vk-s)/k 



(58) 



Hence the explicit solution in the non-ionized region £ > 
is 



<t(x, (>0,t) = = p(x, ( > 0, t), 



(x,C>0,f) = < + <5 



-fcC 



a ikx-\-st 



(59) 
0(5 2 ). 



C. Solution strategy and limits for C, — > — oo 

We aim to calculate the dispersion relation s = s(k) 
for fixed /c. For any s and fe, the solution at £ > is 
given explicitly by (|59|) . This solution determines the 
value of the fields (|58|l at ( — as a unique function of s 
and fc. The expression ()58|l is the initial condition for the 
integration of Ij45(l towards C — > — oo. The requirement 
that the solution approaches a physical limit at ( — > — oo 
has to determine s as a function of k. According to a 
counting argument, this is indeed the case, as will be 
explained now. 

First, the limiting values of the fields at £ = — oo are 
comparatively easy: the total charge vanishes, hence <j\ 
and pi approach the same limiting value <j\ — > er-f and 
pi — > erj~, and the electric field vanishes, hence ipi — > 
and </)i — > 0. Here the limiting values at ( — > — oo again 
were denoted by the upper index ~ as in (|34[) . 

Second, the eigendirections are determined by lineariz- 
ing the equations of motion l)45|l about this asymptotics. 
In a calculation similar to the one from Sect. II. F, one 
derives for £ — ► — oo 



Pi | c^-°° 




+a 3 e 



with the free parameters 01^2,03,04 and a 1 and the 
eigenvalues 



A 1 = ^^ = A + A 2 , A 2 = - (61) 

v v 



For positive s and fe, all eigenvalues Ai, A 2 and fc are 
positive except for the fourth one —k. Hence the first 
three eigendirections approach the appropriate limit for 
£ — > — oo, while the fourth one does not. Therefore a 
solution can only be constructed for 



CI4 = 0. 



(62) 



and A from Eq. 1331) . 



This condition determines the dispersion relation s = 
s(fc) when a solution of ll4"5|) and is integrated to- 
wards £ — > — oo. 



IV. CALCULATION OF THE DISPERSION 
RELATION 

Having set the stage, the dispersion relation is now 
first evaluated numerically for = — 1. Besides an 
expected result for small k, this investigation has deliv- 
ered a previously unexpected result for large k. Based 
on these numerical results for fixed Eoo, we were able 
to derive analytical asymptotes for small or large k and 
for arbitrary < 0. We also understood the physical 
mechanism driving this asymptotic behavior. The sec- 
tion contains the derivation of our numerical results and 
of our analytical asymptotes and their physical interpre- 
tation. 



A. Numerical results for arbitrary k and = — 1 

The problem is to integrate the equations for the 
transversal perturbation (|45|l for fixed k and guessed s 
from the initial condition H58|) at C = towards decreas- 
ing £. In general, the boundary condition l|60|l with 162|) 
will not be met, so s has to be iterated until 04 « 0. 
When the condition is met, the solution does not diverge 
for large negative Cj otherwise it does. When passing 
through the appropriate s = s(fc), the sign of the diver- 
gence changes. This is how the data of Fig. 2 with their 
error bars were derived. 

For the numerical integration, the ODEPACK collec- 
tion of subroutines for solving initial value problems was 
used [2^ to solve the seven ordinary differential equations 
for the unperturbed problem (|1U|I - H12[) and the perturba- 
tion H45fl - H46[) simultaneously. The unperturbed solution 
has to be calculated since it enters the matrix l|46(l . 

However, the numerics can not directly be applied 
to the problem in the form l|45l) (|46|l because the ma- 
trix contains apparently diverging terms proportional to 
l/(v + E(()) for C — > 0. Therefore the behavior of the 
solution for £ — > has to be evaluated in a similar way 
as in Sect. II. E. With the ansatz 

a 1 (C)=a 1 (0-) + C 1 C + O(C 2 ) , 

Pl (C) = Mo-) + c 2 c + o(c 2 ) , 
MO = ^i(o-) + c 3 c + o(c 2 ) , 
MO = <Mo-) + c 4 c + o(c 2 ) , (63) 
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where cti(0 ) etc. are given by Q58fl. the parameters Ci 
become 



C-) = -sa 



C, = a -k 



ff 
s + f 
ff 



s + f 



+ / + /'] , 
, C± = s, (64) 



C 2 + (a + W/2) C 3 + s(ua/" + va'f/2) 
2 + s/f 



In the numerical procedure, the explicit solutions l|28[l - 
ignjl and are used until C = 10~ 5 , then the 

differential equations are evaluated. 
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FIG. 2: Dispersion curve for Boo = —1, hence v = 1. The big 
figure shows the numerical data with error bars and the two 
analytical asymptotes for small and large k (lines). The inset 
shows the same data (squares) in double-logarithmic scale 
with the same two analytical asymptotes. 



The numerical results for the dispersion relation in a 
field Eoo = — 1, i.e., for a shock front with velocity v = 1 
are shown in Fig. [21 It can be seen that the dispersion 
curve for small k grows linearly, but then turns over and 
finally for large k saturates at a constant value. 



B. Asymptotics for small k and arbitrary Eoo < 

We first derive the asymptotic behavior for small k for 
an arbitrary far field Eoa < 0. 

When the equations of motion l|45|l and (|46|l are eval- 
uated up to first order in k, <j>\ decouples, and we get 



01 \ 








= N S - 


r 









sd ( a/(v + E) 
sd^p/v 




-0(fc 2 ) , 
(65) 



where 



N S (C) = 



/ s + 2a - f - p -o d ( a - af \ 



E 



V 



1 



E 



-1 



+ E 
-of 



'O(k') 



/ 

(66) 

is the truncated matrix M Si fc(£) l|46|). The matrix N s for 
s = reduces to the matrix No from Eq. (|37f) — this 
fact will be instrumental below. The fourth decoupled 
equation reads 



dc_4>i = tpi 

The boundary condition l|58|) reduces to 

fsf/(s + f) 
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0(k 2 ) 



and 



l(0) 



vk 



(67) 



(68) 



(69) 



Now compare the mode (cte, Pe,iPe) of infinitesimal 
change of far field E^ from Eqs. l(4TJ|) and l(4*T|l to 

the present perturbation mode in the limit of small k. 
After identifying 



(<Ti,pi,ipi) = (sa-E,spE,sip E ), 



(70) 



the equations and boundary conditions for the modes 
are identical in leading order of the small parameter s. 
Therefore the two modes have to become identical in 
the limit s f(v) < v. Integration over ipE yields for 

the electric potential 0b(O) — 4>e(— oo) = /_ dx i/je(x). 
This expression has to be of order unity since all other 
quantities are of order unity. But this implies that <j>\ (0) 
due to (|7UI) has to be of order s. Now compare the result 
for 0i (0) in ^ which appears to depend in a singular 
way like 1/k on the small parameter k. But for small k 
and s the expression (vk — s) / k indeed can be of order s, 
namely if 



= vk + 0{k 3 ) for k < a(v) 



(71) 



This fixes the dispersion relation s = s(k) in the limit of 
small k. The asymptote (|71f) is included as a solid line 
in Fig. 2. 



C. Physical interpretation of the small k asymptote 

This result has an immediate physical interpretation: 
for small k, the wave length of the transversal perturba- 
tion 2-w/k is the largest length scale of the problem. It 
is much larger than the thickness of the screening charge 
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layer that is shown in Fig. 1. Therefore on the scale 1/fc, 
the charged front layer is very thin and has the charac- 
ter of a surface charge rather than of a volume charge. 
This surface is equipotential according to l|59|l in linear 
approximation in the perturbation S, since 



(x,( = 0,t) = 6 



vk 



ikx-\-st 



k 



0{8 2 ) 



0{Sk) + 0(S 2 ), 



(72) 



if we insert the dispersion relation s — vk from (|T 1 1) . The 
corresponding electric field ahead of the interface is 

E[x,( = 0+t) = -(v + 5vk e lkx+st ) ( + 0(5 2 ) (73) 

in the same approximation. The small k limit of the 
ionization front therefore is equivalent to an equipotential 
interface at position C — 0, i.e., at a position 



z(x, i) =vt + 8 e 



ikx+st 



(74) 



in the rest frame z (14211 . Its velocity in the z direction is 
therefore 



v(x, t) = dtz(x, t) 



Svk e lkx+s \ 



(75) 



where s = vk was inserted. Comparison of l|73|) and l|75|) 
shows that the interface moves precisely with the electron 
drift velocity v — —E within the local field E. 

We conclude that a linear perturbation of the ioniza- 
tion front whose wave length is much larger than all other 
lengths, has the same evolution as an equipotential inter- 
face ((/> = const.) whose velocity is the local electron drift 
velocity v = V</>. It exhibits the familiar Laplacian inter- 
facial instability s oc k. 



D. Asymptotics for large k and arbitrary E x < 

For large wave-vector k, the numerical results for the 
dispersion relation s(k) in a field E^ = — 1 approach a 
positive saturation value. We will now argue that the sat- 
uration value is given by s(k) = f(E oa )/2. This asymp- 
totic value, which for \Eoc\ = 1 equals e _1 /2 = 0.184, is 
included as a solid asymptotic line in Fig. 2. 

When the electron and ion densities remain bounded, 
the equations with the most rapid variation in (|45[) - H46|1 
for k 3> 1 are given by 



<9 c </>i 



k 2 ^ 



k 2 E(C)+O(k ), 



= fa 



(76) 



On the short length scale 2-7r/fc, the unperturbed electric 
field for C < can be approximated as in i|30[) by 



E(0 = -v-f(v)C + O(( 2 ), 



so the equation for <pi becomes 
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d^ 1 =k 2 {<p 1 -v-f(v)C). 



(77) 



(78) 



The boundary condition i|58[l fixes </>i(0) = (vk — s)/k 
and "01(0) = <9(^0i = s. The unique solution of l|78l) with 
these initial conditions is 



0i (0 



f(v)C 



m 

2k 



M 



f(v) - 2s 
2k 



(79) 



for C < 0. Now the mode e would increase rapidly 
towards decreasing £, create diverging electric fields in 
the ionized region and could not be balanced by any other 
terms in the equations. Therefore it has to be absent. 
The demand that its coefficient (f(v) — 2s) /2k vanishes, 
fixes the dispersion relation 



m = 



/(«) 



for k > a(v) , (80) 



which convincingly fits the numerical results for large k 
in Fig. 2. 



E. Physical interpretation of the large k asymptote 

Also for this result a physical interpretation can be 
given. First note that the z component of the electric 
field on the discontinuity is 

E z (x, C = 0, t) = -(« + 6 s e lkx+st + 0(S 2 )) (81) 

with s = f(v)/2. This is easily determined from cither 
(|59"|l or (03) ■ Reasoning as in l(73* )l -l(7H |l . we again find 
that the shock line of the electron density moves with 
the local electron drift velocity — as it should. 

Second, one needs to understand why the electric field 
on the shock line takes the particular form 181fl . In the 
frame £ = z—vt of the unperturbed front l|42|) . the electric 
field at the discontinuity is 



E z (x,£ = A,t)=- 



A + 0(A 2 



(82) 



from 



where its position deviates with A(x,t) = S e lkx ^ 
the planar front. 

In linear perturbation theory, the amplitude S of the 
perturbation has to be much smaller than its wave length 
27r/fc. Since this wave length 27r//c now is much smaller 
than the width of the front, the linear perturbation A 
explores only a small region around the position of the 
shock front. In this region, the electric field of the un- 
perturbed front is according to (|30f) approximated by 



f(v)t; + o(e)) 



for 
for 



(83) 



Therefore the electric field (j82(l is just the average over 
the behavior ijSSjl for £ > and £ < 0. This spatial 
averaging is enforced by the harmonic analysis of linear 
perturbation theory that will suppress different growth 
rates of positive or negative half-waves of the perturba- 
tion. 
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F. A conjecture for the large k asymptote 

We therefore conjecture: if the electric field of an un- 
perturbed front is 



#o(0 = 



of- 



-o(e)) 



for 
for 



£<0 



(84) 



near the position of the discontinuity £ = 0, then a linear 
perturbation of this discontinuity with large k will grow 
with rate 



(85) 



If true, this behavior would have a stabilizing effect 
on large k perturbations with growing curvature of the 
fronts, since the electric field decays in the non-ionized 
region ahead of a curved front, therefore b < 0. 



V. CONCLUSIONS AND OUTLOOK 

We have studied the (in)stability of planar negative 
ionization fronts against linear perturbations. Such per- 
turbations can be decomposed into transversal Fourier 
modes. We have determined the dispersion relation 
s = s(k) shown in Fig. 2 numerically for a fixed field 
= — 1 far ahead of the front, and we have derived 
the analytical asymptotes 



\Eoo\k £arfe<o(|£oo|)/2 
l-Bool afl-EooD/2 forfc>a(| J E oc |)/2 



(86) 



for arbitrary E&0 < 0. Since we have studied the minimal 
model, there is only one inherent length scale, namely 
the thickness of the charged layer as shown in Fig. 1. 
This thickness is approximated by l/a(-B c>0 ). The wave 
length 1/fc of the Fourier perturbation therefore has to be 
compared with this single intrinsic length scale l/a(£ , oc ) 
of the problem. 

A specific property of our calculation is the expansion 
about a discontinuity of the electron density. Therefore 
we work in a coordinate system £ = z — vt — fi e lkx + st 
(|42|1 that precisely follows the position of the discontinu- 
ity, and we explicitly distinguish in all calculations the 
non-ionized region £ > from the ionized region £ < 0. 
For the non-ionized region £ > 0, there is an exact an- 
alytical solution l|59|) for any s and k which determines 



the values of the fields at £ = as given in lj55|) . Eq. lj5"5|) 
serves as an initial condition for the integration towards 
£ < 0. The approach towards £ ^ — oo according to l|rjU|) 
and I|tj2f) determines the growth rate s as a function of 
k. In general, this calculation has to be performed nu- 
merically with results as shown in Fig. 2. The limits of 
small and large k can be derived analytically. For small 
k, we can identify the perturbation mode with the mode 
of infinitesimal change of E^. For large k, the growth 
rate corresponds to the evolution of the discontinuity in 
the unperturbed electric field averaged across the dis- 
continuity. Both limits therefore have a simple physical 
interpretation. 

The aim of the work was to identify a reg ularization for 
the interfacial model as suggested in jaUj and treated in 
[Iol| . Indeed, we have found that a Fourier mode for large 
A: in a far field Eqo — —v does not continue to increase 
with rate s = vk, but saturates at a value s = f(v)/2. 
Still this is a positive value, and whether this suffices 
to regularize the moving boundary problem, remains an 
open question. 

Besides this one, future work will have to investigate 
two more questions. First of all, there is the "simple" 
possibility to extend the model by diffusion. Diffusion 
is certainly going to suppress the growth rate of Fourier 
modes with large k as our preliminary numerical work 
indicates. But there is also a second more subtle and 
interesting possibility: the growth rate of Fourier pertur- 
bations with large k could change for a curved front, as 
we have conjectured in Section IV. F. There we have ar- 
gued that the saturating growth rate s = f(v)/2 results 
from the average over the slope —f(v) of the field in the 
ionized region and the slope of the field in the non- 
ionized region. For a curved front, the electric field in 
the non-ionized region will have a slope of opposite sign 
that is proportional to the local curvature. We there- 
fore expect the growth rate of a perturbation to decrease 
with growing curvature. These questions require future 
investigation. 
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